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In this paper we further explore and develop the quantum continuum mechanics 
(CM) of [Tao et al, PRL103, 086401] with the aim of making it simpler to use in 
^ . practice. Our simplifications relate to the non-interacting part of the CM equations, 

r-| ■ and primarily refer to practical implementations in which the groundstate stress ten- 



sor is approximated by its Kohn-Sham version. We use the simplified approach to 
directly prove the exactness of CM for one-electron systems via an orthonormal for- 
mulation. This proof sheds light on certain physical considerations contained in the 
CM theory and their implication on CM-based approximations. The one-electron 
proof then motivates an approximation to the CM (exact under certain conditions) 
expanded on the wavefunctions of the Kohn-Sham (KS) equations. Particular atten- 
tion is paid to the relationships between transitions from occupied to unoccupied KS 
orbitals and their approximations under the CM. We also demonstrate the simplified 
CM semi-analytically on an example system. 
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I. INTRODUCTION 



The study of the continuum mechanics (fluid dynamics) of quantum electron fluids is al- 
most as old as quantum mechnics itself. The Thomas-Fermi model and Madelung dynamics^ 
are two very early examples. Unfortunately standard approaches to Continuum Mechanics 
become too inaccurate or too complex when applied to larger, real systems such as molecules 
and atoms, although a recent development^ makes atomic systems more tractable. Recent 
work^"''' on Continuum Mechanics (CM) in a moving Lagrangian frame has led to the de- 
velopment of a sophisticated approach for describing linear perturbations to many-electron 
systems from groundstate properties only. This approach appears able to bridge the gap 
between speed and accuracy required by modern ab initio calculations. 

The CM provides an efficient^ alternative to full time-dependent density functional theory 
(tdDFT)^"^^ calculations. In its general form it can, in principle, be used to evaluate the 
transition frequencies and currents of a real many-electron quantum system using the one- 
and two-particle density matrices. In the formalism presented here we restrict to a more 
limited form that takes, as input, groundstate properties obtained from a Kohn-Sham® 
calculation, and approximates small changes to the groundstate via a continuum approach. 
This restriction to groundstate KS properties only comes at the expense of having to deal 
with higher order mixed derivatives (up to four derivatives with three indices). Initial 
indications suggest, however, that it is both tractable and vahd in both model systems, ^'^ as 
well as the difficult and geometrically very different case of two interacting two-dimensional 
electron gas layers.'^ 

Like Madelung dynamics, the recent CM approach describes the behavior of the fluid 
displacement vector u from which the current and changes to the density can be described. 
As the independent-electron density response Xo of a system can be obtained from w, a "CM- 
dRPA" correlation energy functional^ has been developed that uses the direct Random Phase 
Approximation (dRPA) but bypasses the need for unoccupied orbitals by working in the 
CM directly via u. The CM-dRPA can be considered a "third-rung" functional according 
to the "Jacob's Ladder" classification scheme of Perdew et al}^ The functional involves 
calculating the bare response via the CM scheme and solving for the interacting response xx 
under the dRPA where additional interactions are treated at time-dependent Hartee level 
via Xa = Xo + ^Xo^Xx- The response functions Xo and xx can then be used to calculate 
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correlation energies. Here A is the strength of the interactions and must be integrated over 
to obtain the kinetic contribution to the correlation energy. 

One particular area where it is hoped that the CM will prove broadly useful is in the 
evaluation of van der Waals physics, where long-range correlation is important. Here local- 
density based techniques^^"^'^ and even vdW-adapted approximations like the vdW-DF^®~^^ 
run into difficulties (see ref. 7 and ref. 22 for discussion) because of their use of local or 
pairwise approximations. 

The strong theoretical relationships between the CM and the KS-like system it approx- 
imates, arising from its derivation from a formal, moving Lagrange frame,^ is somewhat 
hidden in the current prescription, especially to those used to working with orbitals and the 
Schrodinger equation. While relationships can be established (such as sum rules) between 
u and common groundstate properties of interest, they do not always come naturally from 
the original formulation. 

Via adaptation and exploration of its theoretical form the CM offers a wider scope for 
further investigation, as follows: 

Firstly, in their original form^ the CM equations were complicated. In the first application 
to vdW physics,^ a more compact, simpler and more symmetric but equivalent form was 
used (equations 3-5 of ref. 7) . The derivation of this simplified form was not given in ref. 7, 
and is presented here for the first time along with a discussion of the particular version of 
the groundstate stress tensor that makes this simplified version of the CM theory possible. 

Secondly, the standard displacement vector u has some undesirable properties in Coulomb- 
localised systems such as atoms and molecules, including divergence in the tail of the density 
distribution. To deal with this we reformulate the CM in terms of an "orthonormal dis- 
placement vector" ^ = VnPu (touched on in equation (68) of ref. 6) where nP is the electron 
density of the groundstate. This alternate approach opens the method up to expansion in 
basis sets with decaying tails, such as the widely used^^ Gaussian- Type Orbitals (GTOs) and 
Slater-Type Orbitals (STOs) which are only vahd in the expansion of decaying functions. 

Thirdly, this reformulated CM can be used to demonstrate, directly from the Schrodinger 
equation, that the CM gives the exact bare, linear response xo for one-electron systems. 
While this relationship is established and demonstrated in earlier works (see equations (41)- 
(46) in ref. 4, equation (98) and Appendix D in ref. 5) a proof direct from the Schrodinger 
equation has not so far appeared. This relationship is very important, having as a con- 
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sequence that properties dependent on xq such as the (dRPA) asymptotic van der Waals 
interaction between atoms, are exactly reproduced by CM-dRPA theory for one-electron 
systems (and two-electron systems with equal groundstate densities of spin up and down 
electrons). The direct proof demonstrates how the CM relates exactly to the first-order 
change of the one-electron orbital/wavefunction. 

Fourthly, the one-electron case motivates a second, approximate reformulation of the 
tensor CM into a scalar system. Here the displacement u is approximated as the gradient of 
a scalar function s. This approach simplifies the CM equations at the expense of accuracy 
in general systems but is exact for one-electron and one-dimensional systems. By then 
expanding VnPs on the set of KS orbitals we also uncover some of the physics of the CM 
in the asymptotic tail regions of the density, where the electrons behave like a one-electron 
system. 

Finally we illustrate the work on an example system: a one-dimensional, non-interacting 
Harmonic potential model. Here the one-electron case can be solved analytically while many 
of the terms in the many-electron case can be solved for analytically or using near exact 
quadrature, minimising numerical error. The scalar approximation is an exact reformulation 
of the CM in this example. 

A. Notation 

In this paper we work entirely in atomic coordinates where rUe- = h = e^/(47reo) = 1 
such that energies are in Hartee and distances in Bohr radii. We treat energies as frequencies 
with the division by h implicit. 

Greek sub/super-scripts are used to refer to Cartesian {x, y and z) coordinates and are 
summed over if repeated. We use the derivative operator notation [9^/(t*)] = [9/(r)/9r^]. 

Cartesian tensors are written with sans serif letters (eg. T), while Cartesian vectors 
appear in bold (eg. v). Their elements are typically given as T^,^ and respectively. More 
general matrices use double-line Roman letters (eg. M) and should be considered square 
unless otherwise noted. 

Operators involve a hat eg. O. If a derivative appears in an operator it can be considered 
to act entirely to the right unless surrounded by square brackets, but will act through other 
brackets. Thus [daVa] = [V • v], dj = [d^f] + fd^ and (5^/)c/ = [d^,f]g + f[d^g\ + fgd^. 
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Comma-led subscripts will sometimes be used to represent derivatives [d/^fv] = fu,ij,- 

In the context of KS orbitals j or k can typically be any orbital but we reserve i for 
occupied orbitals only and a for unoccupied orbitals only such that J2i = J2i occ Sai = 

T- y 

z__/i occ ^-^a unocc 

II. ORIGINAL VS OPERATOR FORMS OF CM 

Let us define a Kohn-Sham (KS) system with a KS potential V^^ (approximate or oth- 
erwise). The one-electron Hamiltonian is thus 

|_lv^ + y^^(r)}v',(r)=6,V',(r) (1) 

where tljj{r) — {r\j) is a one-electron orbital which are orthonormal under J drtlj*{r)ipk{i') = 
{j\k) — 5jk and where is its Kohn-Sham eigenvalue. 

The groundstate density n° is, in general, defined via n(r) = p{r, r) where the one- 
electron density matrix p(r, r') is discussed later. In a KS system the density becomes 

n'{r)=J2f^rj{r)Mr) (2) 
j 

where fj is the occupation of orbital \j) defined as 1 for orbitals with ej < ep and otherwise 
where is the Fermi energy, ep should be chosen to ensure that J n^{r)dr — where A^e 
is the total number of electrons. 

A. Kinetic stress tensor 

For the stress tensor T corresponding to a onc-clcctron density matrix p(r, r') in the 
absence of a current we choose the following definition, which we motivate and derive in 
Appendix A: 

TAr) =\{dA + ^^d^)pir,r')\r=r' 

- ^d^d^nQ{r). (3) 

This particular definition seems to allow for the most compact presentation of the CM 
eigen-equations (18). 

5 



Inserting the density matrix from the independent-electron Kohn-Sham groundstate we 
obtain 

t;. h - (4) 

^^^E/^([^mV';][5.V',] - rAdM) (5) 

where terms are functions of r only. Unless otherwise noted we subsequently restrict our- 
selves to this form. The kinetic stress tensor T° is real symmetric [T'°j^(r) = T^^{r)\ and 
obeys the groundstate force balance condition [9aT°^] = — 

In Appendix A we discuss the origins of the form used here and reasons for the difference 

between our forms of T and T" and the form defined in equation 17 of ref. 5 [equation 13 
here, leading to our equation 14 for the Kohn-Sham groundstate.] 



B. Linear CM 

In general, there is a linear relationship between the first-order time-dependent fiuid 
displacement vector u(r,i), and a small change to the external potential defined via its 
force density per electron F^^^^{r, t) = {r)W^^'^^ {r , t). Here u is related to the electron 
current through^^ 

j(r-,t)=9tn°(r)u(r,i), (6a) 
dtn\r,t) = -V ■j{r,t) (6b) 

from which it follows that a solution is 

n\r) ^-V ■n%r)u{r). (7) 

Linear CM under a Kohn-Sham groundstate provides this relationship approximately in the 
following form for a time-periodic potential of frequency uj^'^ 

Mu + F^^^" + F^"^"* + F^^"* , (8) 

where F^^"* includes the internal dynamic electron-electron interactions and is discussed 
later. The remaining terms are: F^^™ — — Kw where^ 
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is an Hermitian operator, a full derivation of which appears in appendix B; And F ° — 
n°[V (g) Vy^^] • u. Putting them together in the form F^^°* + F^^'"" = Ru leads to 

- R)u =F^^"* + F""\ (10) 
R,^ =^°K?' - (11) 

where all terms but ou vary with r. The operator R is manifestly Hermitian, and can be 
shown to be positive definite. 

The original papers^'® establishing the linearised CM theory give F^^^^ [our (9)] in a less 
compact but exactly equivalent manner. The expression for F^^^^ appearing in equation 14 
of ref. 5 (equation 53 of ref. 6) is (where we have rearranged the order) 

-F^i^'° =a„(2f;,t/,, + f °„a^«,) 

+ { [a.nV/. + [dy% - df^n^d,} V • u 

+ {[W]C/^, - d^{[dan']U,a)} . (12) 

Here t/^jy = ^[d^u^j, + dfj^u^,] and their kinetic stress tensor is [equation (17) in ref. 6] 

-\s,,[V'no{r)]. (13) 
For an independent-electron Kohn-Sham groundstate this takes the form 

= ^J2Md,rjWj] - \s,A^'n'] (14) 
j 

where terms depend on r only. The second components in (13) and (14) differ from ours 
[compare (3) and (4)]. Although it is not immediately obvious, the force (12) is equivalent 
to an Hermitian operator acting on u. By contrast the Hermiticity of our operator (9) is 
apparent on inspection. 

Equation (9) defines the same force as (12) but in a much simplified manner. Much of 
this simplification comes from the different choice of kinetic stress tensor, given by (3) and 
(4) in our work and by (13) and (14) in refs 5 and 6. In fact there is some ambiguity in the 
definition of the quantal stress tensor. In Appendix A we use the Wigner distribution to 
motivate our present form, which is then used to derive (9) in Appendix B. 
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III. CM RESPONSE FUNCTIONS 



In (8) takes into account the small change in the Kohn-Sham potential arising from 
the change in density, while Ru deals with kinetic and potential physics. The dRPA is 
equivalent to setting 

F"°*(r) =[Qw](r) (15) 

/Hr' 
^^[V'-n°(rOw(r')] 

= J dr'Q{r,r') ■ u{r') (16) 



where 



Q,.{r,r') =n'{r)n\r'' 



'r — r 



(17) 



In the absence of an external field F^^^^^ = we can find eigen-mode pairs Vl^ and Un 
(or VLnx and Unx) through solutions of 

Q%n°UN =Rwiv, (18) 
n%^n°UNx ^Runx + AQwjvA (19) 

where (19) includes the internal interactions AQ at coupling strength A while we use the 
short-hand Qat = and Un — Unq for the non-interacting case. Because R and Q are 
Hermitian and is symmetric and positive definite, eigensolutions can be found that obey 
the orthonormality condition 



drn u^;^ ■ umx = Snm- (20) 

The set {unx} is also guaranteed to be complete over a finite basis when R-I-AQ is represented 
in the same finite basis. Furthermore the eigen- values can be shown^ to be positive. 
Typically we also sort the modes such that ^at+ia > f^iVA where > 1 labels the CM mode 
with the A^th lowest energy. The displacement Unx corresponds to a transition density mode 
(the meaning will become clearer later) defined as 

dNx{r)^-V ■n''{r)uNx{r). (21) 
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If the external force density ir^^xt reintroduced we can expand the solution of (8) at 
interaction strength A on the basis {unx} such that u — cnxUnx- Here 

^n;;^^ ^^'^ 

when F^^^^ is time-periodic with frequency u. The change in density (7) can thus be 
expanded on (21) as 

n\{r,uj) ^^CNxdN\{r), (23) 

AT 

where the sum is over all eigen-solutions. 

The density response x\ of ^ system is defined as the change in density in response to a 
S{r—r') potential at a frequency uu with internal interactions at strength A. This corresponds 
to an external force F^^""^ = n^V5{r-r')e-''^^ and internal force F^^"* = AQw(a;)e-^'^*. Thus 
the response takes the form Xxit) = XA(i^)e~"^* where xa(^, i^) = — n ~^^^w^~^^^~^ ■ 
Typically it is easier to work with responses at imaginary frequency uj = ia such that 

where (24) uses the solutions of (18) to calculate the bare (A = 0) response while (25) uses 
(19) to solve directly for the interacting response. The 3? is unnecessary as the sum itself 
can be guaranteed real but may prove useful in some situations. 

While xx defined by (25) has useful formal properties its direct evaluation may be nu- 
merically difficult and can be avoided. Unless otherwise noted we henceforth set F^^"^^ — 
and deal with internal interactions (when required) in a less direct, but precisely equivalent 
manner (as discussed tangentially in ref. 7 and in Appendix C of this manuscript). 



IV. RELATIONSHIPS TO KS RESPONSE 

In a Kohn-Sham system with orbitals i'j{r) = {r\j), the exact bare response takes the 
form^ 

Xf(r,r^.a)=-5R 5:^^:1^ (26) 
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where ^ai is a normalised transition density between unoccupied orbital \a) and occupied 
orbital |i) while is the transition frequency defined by 



7„i(r) ^^/2n^^l{r)'il^i{r), (27a) 
^ai =ea - Q > 0. (27b) 

As noted in Sec. I A, i is summed over occupied orbitals only and a over unoccupied orbitals 
only. 

There is a transition current density associated with \a) and |i) which takes the form 

Jarir) =^[^a(r)V^*(r) - ^*(r)V^„(r)] (28) 



and where iV • j„j = Qai'^i'^l — s/^ai/'^lai- Since {wat} is complete (at least within a given 
finite basis) and orthonormal under (20) we can expand j^^ as 

jai = ^aiNn^UN, K^iN = / ^Tu]^ ■ J^^. (29) 

TV 

Taking the gradient of thus provides the following relationship between the KS transition 
densities and the CM density modes djv 

TV V »v 

so that any can be expanded in {dj^}. Unfortunately the converse cannot be guaranteed 
except in the trivial one-electron case. In molecular systems we can make both cLn and ipj 
real so that 

dN = ^ iNailai + ^ Im'ili'i (31) 
ai i<i' 

where both i' and i are occupied while a is unoccupied. 

Certain exact sum rules [equations (81)-(83) further discussed in appendix E of ref. 6] 
provide some further restrictions on the various coefficients. Since the f- and third-moment 
sum rules are satisfied by the CM it follows^ that 

1 = ^ haiN, ^^AT = ^ haiN^ai (32) 

ai ai 

where Qai are the Kohn-Sham transition frequencies if the CM equations involve the KS 
density n° and stress tensor T*^. Here the non-negative weights haiN are defined via 

iT-aiN — 7^ l^-J-jj 
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where KaiN is defined in (29). More generally iox N ^ M these take the form 

^ 2Kl,^K^,M ^ ^^^^ ^34^^ 

ai 

E ^%^^a. = ^NM^l (34b) 

which come from inserting (30) into (26) and comparing the leading two powers of 
with (24). 

The second sum rule makes the relationship between Qjv and the KS transition frequencies 
clear. We may also consider to be an approximation to collections of the transition 
densities, with errors hopefully minimised by the sum rules and exact properties although 
no direct expansion exists. As discussed later these approximations become exact for one- 
electron (or '\\) systems. 

It is also worth noting that the lowest CM transition frequency Jli can never be less than 
the transition frequency between the highest occupied- and lowest unuccopied- molecular 
orbital ^lh — ^l — ^h- In the non-degenerate case the equality follows if and only if Hlhi — 1 
with all other han zero. To prove the inequality we note that Q^ai ^ ^lh and thus 

^l>Y.'^-i^^lH>^lH- (35) 

ai 

If hLHi < 1 then hLm = 1 - T^ai^LH ^aii and 

= E ha^,{nl^ - nl^) > 

ai^LH 

since han > and fl^^ — Qj^jj > 0. Thus the equality only holds if haa = 1. 



V. CORRELATION ENERGIES MADE SIMPLE 

Once we have obtained and Jljv, the CM approximation to the correlation energy can 
be calculated.^ We define the Coulomb projection matrix W with elements 



W; 



NM 



= j drul,{r)-[QuN]{r) (36) 

ardr'^JiP^ (37) 
\r — r'\ 
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and L with elements Lmn — ^mn^'n- Through the working in Appendix C we can show 
that the correlation energy is 

d(T 



Jo Jo TT 



X -Tr 
2 



W 



(72 + L + AW (72 + L 



(38) 



or we can use the Furche-like^^ integrated form 



N 



N 



NN 



2Q% 



(39) 



where arc the eigenvalues of R + Q or L 



L As discussed in Appendix C the two 



diagonalisations arc formally equivalent but experience in similar techniques suggests that 
working in the transition densities of the bare response will allow for better convergence. 
In practise diagonalising L + W is expected to be faster and numerically more reliable and 
robust. 

Using the eigenvalues of L+ W has a further advantage: we can use a perturbative solution 
to find the eigenvalues of L + W if flj^ 3> Wjy where Wn = J2m I^a^mI (see Appendix C for 
details). We define an N* such that fl% > KWn^N > N* where K is sufficiently large. We 
then solve the reduced N* x N* eigen-equation L* + W* to obtain and calculate 



A- 

Ec^-Y,[n*^-n^{i + M]- E 



N 



(40) 



N=l 



N>N* 



where Pn^^- 



VI. ORTHONORMAL DISPLACEMENT 

If we consider the orthonormality condition (20) of the CM eigenmodes in a bound system 
we can sec that Un may be permitted to grow as |r| — )> 00 provided VnPuj\f decreases. In 
atomic and molecular systems all valid solutions will, in fact, grow exponentially due to 
the asymptotic form of the orbitals. This suggests that an orthonormal fluid displacement 
^ = VnPu will be a more natural quantity to use in these systems as it is guaranteed to 
decrease. Here the orthonormal modes are 

= V^Un (41) 
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where / dr^^ • £m = ^nm 

While u — X^^yCivWAr has a well-defined phsyical meaning (the fiuid displacement), ^ = 
X^cat^jv is somewhat harder to interpret. However some insight can be gained if we define 
the groundstate quasi-orbital = \fnP and its first order change . Since = 2^^°^^ = 
-V • it is clear that 

= + (42) 

and thus ^ is related to the perturbation of the quasi-orbital. It is also related, via (6), to 
the groundstate properties and j through 

j =dt'^% = - V*°^. (43) 

We must calculate through the CM equations (18) 

=R'^'««, (44) 

where R is defined in equation (11). The change in density of a given mode becomes 

dr,{r) = - ^0 (V + r7) • Civ (46) 

where 

is the logarithmic gradient of Equation 46 can be used in (37) to calculate the matrix 
elements of W for use in correlation energy calculations 
Inserting (11) into (45) gives 

+ K' + \^o9^9^^9^9,^^. (48) 



Using the derivative operator identity 



J_a -d — =— (A9) 
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and defining t^,^ — T^^^/nP, allows us to convert (48) into the following succinct reformulation 

- {da + Va)taAdi^ - VlJ,) 

- (d„ + r]^)tan{da - Va) (50) 

which we can use in (44). Here 

S - (V + ry) • (V - r/) 

H^aVa] + VaVa = ([^VV]/n° - TjaVa)- (51) 

All functions appearing in these expressions depend on groundstate orbital wavefunctions 
and their derivatives only. Expanding in terms of occupied orbitals they are 

S - VaVa. (54) 

Here the force balance equation V^^ = —{da + 2ria)taf^ replaces the usual nPV^^ — —daT^^. 
In a one-electron system with occupied orbital ip and energy eo these reduce to 

Va V = -^^VnVi^ - -j^), (55a) 

5=2(\/^S-eo). (55b) 



It is worth noting that, for molecular systems with Coulomb-like nuclear potentials the 
outermost tail is dominated by one-electron-like behaviour and: 

1. The denominators of t^^, rja and 5' are densities and thus everywhere positive, 

2. lim^^oo V(^) ^ 0' 

3. limj._^oo |'7(»')| ~ W~'^^h\ where €h is the KS eigenvalue of the highest occupied 
orbital, 
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where r is the displacement from the center of the highest occupied orbital and we ignore 
leading terms that decay exponentially. In large molecules these expressions may also hold 
true closer to nucleii A with €h replaced by the local e^. 

VII. ONE-ELECTRON SYSTEMS 

Following equations (41) and (45-46) of ref. 4, it is possible to show that the non- 
linearised CM formalism is equivalent to Madelung hydrodynamics^ in a Lagrangian frame 
for one electron systems. If all external fields involve gradients of scalar potentials only, this 
is directly equivalents^ to finding solutions of the one-electron Schrodinger equation (SE). 
Here wc show this equivalence directly from the SE in the linear response limit required by 
the density response Xo- 

We proceed to prove this as follows: i) we derive the relationship between the perturbed 
one-electron wavefunction and its "orthonormal displacement vector" ^^'^ to show that the 
latter is entirely determined by the former (and vice versa up to a trivial phase via the Runge- 
Gross theorem^); ii) we show that finding a free-standing solution of the linear-perturbed 
SE for a onc-clcctron system is equivalent to solving an equation of form cu^^^*^ = R^^ ; 
and iii) we show that R^^ = R as defined in (11). 

In a one-electron system we can set V — V^^ — Eq to eliminate the KS energy of the single 
occupied orbital. Thus the groundstate Hamiltonian takes the form 



where -0 is the only occupied electron wavefunction which we make real. If we apply a small, 
time-dependent external potential V^{t) then we can find a new solution ip'{t) through 



which will be perturbed only slightly from the groundstate solution. We can write the 
perturbed wavefunction via a change to its magnitude and a rotation of its phase such that 




(56) 



[h + V\t)W{t) 



(57) 



(58) 
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or, truncating to first order, '0 + + ii/j(f)^{t). Thus it is sufficient to calculate 

(f)^ and to fully determine the perturbed solution. 

Inserting (58) into (57) and matching real and imaginary components (to linear order) 
gives 

^ljdt(f)\t)^^ljV\t) + h^lj\t), (59) 
M\t) =(V^) • [Vcp\t)] + l^V'cP'it). (60) 

where the second expression relates dtip^ directly to V^^. If we then assume a time-periodic 

external potential V^{t) = V^e''^* it follows that ip^{t) = ip^e'""^ and 0^(t) = and thus 

dt = too. Wc can then use (60) to eliminate ijj'^ and derive the equations governing density 
perturbations in terms of (p^ only. Here 

=li^''Hi^,a + d^ij)d^(p' + tujV\t), (61) 

=2V'V'' = -^ita + da^)da(t>' (62) 
—zu 

where is the linear change in density and wc have used the derivative relationship daip — 
j/jda = ■0,a- By the Runge-Gross theorem^ V' — V + V^{t) is a functional of the density 
n' — nP + ri^it) only in the hnear response regime, and from (62) it is clear that t/j and 
are sufficient to determine electronic properties. 

The Schrodinger current density j of the perturbed, one-electron system is calculated 
through 

f^- =1[^'VV''* - ip'*W] ~ -V'^V(/)^ (63) 

where we use the first-order perturbation expression ip' ^ ip + ip^ + iil^p^ to derive the 
second identity. In generaP^ the displacement u is related to the current via (6) and thus 
iuju^^~ — —Vcf)^ for the one-electron system. 

It then follows trivially that the normahsed displacement = VnPu^^ = ijju^'^ is 
related to 0^ via 

^le- = JLv01 (64) 

—tu 

and, from the Runge-Gross theorem,^ that 0^ can be obtained from . Using (62) the 
density perturbation takes the expected form (43) n} — —V • V'^^^ ■ This completes the first 
stage of the proof. 
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By taking the gradient of (61) and using = {—iu/ip)^^^ [from (64)] we can explicitly 
solve the free-standing (V^ — 0) equation for via 

"^'^ =^5mV'-^/^(V',. + 9.^)^ (65) 
^X" -^\{d, - r,,){y' - S){d^ + r,^)er (66) 

where we used the one-electron specific relationships (55) r]a = '4>,a/'4> ^'^'^ 2h — S — V'^ for 
the second expression. Thus we can define a linear operator 

- - S){d, + V.) (67) 

such that oj'^i]^ = R]^^, Cl^ which is of the same form as (8), as desired for the second 
stage of the proof. 

Finally it remains to be shown that R = R^^ . Following the working in Appendix D we 
show that 

- {da + r]a)t^u{da - Va) 

- {da + r]a)tau{d^ - r]^) 

- {du + r]i^)tan{da - Va) 

^Ri'J (68) 

and the proof is complete. 

Thus we have shown that i) ^^'^ is bijectively (up to a phase) related to the first-order 
solution of the perturbed SE; ii) It obeys cu'^^^'^ = R^'^ for free-standing modes and 
iii) R^^ — R. Thus the governing equation is identical in both cases and it follows that 
= ^. This confirms that the solutions of the CM equations are directly equivalent to 
the solutions of the perturbed Schrodinger equation in one-electron systems, and that the 
CM is thus exact. 
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A. Transition modes of one-electron systems 



In a one-electron system, the density response Xq calculated via orbital transitions (26) 
must be the same as that caclulated via the CM transitions(24). Thus 

lao{rhao{r') ^ dN{r)dN{r') 

^0 = 2^ ^2 ,0 =2^ ^2 ,0 ■ (69) 



where a is summed over the unoccopied orbitals, while ilao = ^a — ^o and 7^0 = ±v2n^^-(/'a. 
Since the equality must be true for all a it follows that the individual numerators and 
denominators of the sums must be paired and we can choose an such that ^jv^ = f^oo 
and rfiv„(r) =7ao(r')- 

Equating (46) and (27) gives 



dN^ = - V'(V + r?) • = V^^^a- (70) 
This has a (non-unique) solution 

U=^^i^-ri)i^a (71) 

which must therefore be a vahd solution of the CM equations. The constant pre-factor comes 
from — (V -I- T7) ■ (V — 77) =2{h — eo) in a one-electron system. 

Tao et al^ tested the onc-clcctron exactness on general s-transitions in the Hydrogen 
atom. We illustrate (70) and (71) on the Is to 2p transition. Here the Is orbital is occupied 
with orbital V'ls = I and energy — e\s — —\ while the l-p^ orbital has ■02p^ = 



ze 1"^ I \J 'Nl'K and e^p^ ~ ~|- '^^^ density is thus — e ^^/tt and 77 = —r. 

Using Vze-''/^ = (^ - \zr)e-''l'^ we can test (70) and (71). With fisp.-i. = | we find 

^2p.-i.=^e-^/2(^f + 2z), (72) 
rf2,.-i. = -^(V-f)-(zf + 2z)e-'-/^ 

1^^18-02^,. (73) 

in agreement with (70) and where 3/4 = 2£l]_s-2%,z as expected. Note that the usual displace- 
ment vector W2p^-is — ^^(-^^ + ^z) grows exponentially with r. 
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VIII. SCALAR APPROXIMATION TO CM 



In many-electron systems, equation (71) is no longer guaranteed to be true. However, as 
we can write any well-behaved vector function as a gradient plus a curl we may set 

Ijv(r-) -*°{Vs^(r) - V X VN{r)} (74) 

where SAr(r) and v^^r) are arbitrary scalar/vector functions with appropriate asymptotes 
(we choose our gauge to ensure V ■ = 0). The form of (71) suggests that we might 
approximate our eigen-solutions by setting 

~^Nir) «^{V - n{r)}Mr) = (75) 

which is equivalent to making the approximation sjv ~ (Pn/"^^ and Vat ~ 0. This approx- 
imation will be exact for one-electron systems and one-dimensional (ID) systems^^ but is 
not true in general. 

The regular CM eigen-equation (44) is equivalent to finding stationary solutions 68/ 6^ = 
of S[$] = I J dr^* • R ■ ^ for vectors ^ satisfying J dr^* ■ ^ = 1. We called these solutions 
and they can be found through fllf^j^ — R£jv where R is defined in (11). Under the 
scalar approximation we restrict our solutions to vectors £ expressible as (V — t/)^ which 
form a connected subspace of all possible ^. We thus look for solutions ^Ar(0) of E[^] which 
are stationary under variation of (p ie. 5£^[^(0)]/50 = subject to / dr^(0) ■ ^(0) = 1. 

The restricted solutions can be found directly by setting S[(f)] = J drcj)* (V + r]) -R- (V — 
ri)4> and setting the constraint to — J dr0*(V + ^) • (V — ri)4> = 1. The general Hermitian 
eigen-equation for (p^ thus becomes 

CiInUn =RUn (76) 

with orthogonal solutions normahsed under J drcffj^N'^cpM — 6nm- Here 

TV"^ = - (V + 7?) • (V - 77) = 5 - (77) 
R^^- (d^ + r]^)R^,{d, - 77,). (78) 

It is obvious that Vl\ > flf since the subspace minimum of S must be equal to or higher 
than its true minimum. 

Such an approximation loses some accuracy and some nice properties of the true CM but 
reduces the problem from a tensor to a scalar. Its exactness in a variety of systems including 
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the one-electron case in any number of dimensions suggests that it might be appropriate for 
vdW calculations in molecular systems. Here the vdW physics are often dominated by the 
asymptotic regions which show one-electron- like behaviour (see Section VIII C for further 
details). 



A. Scalar approximation in KS orbitals 

We can expand (pN — '^jPNjipj in the KS orbitals (or any other complete and orthonormal 
basis set) so that, from (75), 

lAr) = J^P^^i^ - ^i'^)}Mr) (79) 
j 

where ipj are KS orbitals (occupied or otherwise) . In reality we must truncate to the lowest 
iVeas orbitals. The transition density modes are thus (remembering that y/nP = \E''^) 

Mr) = - *° J] PNj{V + ry) • (V - 77) V'^- (80a) 

i 

^^'Y.P^ji^-^'^'i'j (80b) 

i 

J^PNji^i^k - V^^') + S]^^JJ. (80c) 
i 

where the first two properties are true for any basis set but the third property is only true 
for the KS orbitals. 

Projection into the KS orbitals allows some insight into the physical meaning of this 
approximation to be obtained by considering the quasi-orbital Here we can set — 
(2i/\&°) '^jPNjjj where the quasi-transition current 

3 J (V - V) i^j] = ^[*°VV^, - ^,V*°] (81) 

has a similar form to a transition current = ^{ipiS/ipa — i^ayi^i) with the occupied orbital 
replaced by the quasi-orbital of the total density. 

If we pre-multiply (44) by and integrate we find J dr^j^-^j^ = J dr^j^R^^ subject 
to J dr^j^ ■ — ^MN- Using (79) the two integrals become 

/ drl*^-lN^^P*MjPNkNjk: (82) 

jk 

drl^RlN =^P*MjPNkRjk- (83) 

jk 



20 



where 

Njk = j ^r[{d^ - nM*Wa - ila)il^ki (84) 
Rjk = j dr[(a^ - r]^)ij*]R^,[{d, - r;,)^^]. (85) 
Taking derivatives with respect to puj, the eigen-equation (44) thus becomes 

^\^3kPNk —RjkPNk (86) 

subject to the orthogonailty condition Ylijk-^jkPMjPNk — ^nm- Following the details of 
Appendix E we find 

Rjk — 

+ dr[babi,b,rj][bab^b.iJk] 

+ 3 j drib^b^rjKADabM (87) 

where Da = da-rja- 

B. Matrix form of the scalar approximation 

Since all terms involve repeated use of the operator {da — rja) we may simplify things 

somewhat by adopting a matrix notation to represent this operator. For a finite basis set of 
size A^Bas we define the following N^as x -ABas matrices: D^, Tj^^, and W^i, with elements 

Dajk = J drtlj*{da - Va)'4'k, (88) 

Tf,.jk = J dr7P*ti,,7/jk, (89) 

V.-fe = / drrjV^^i^k. (90) 

We can also integrate by parts to obtain V^j/jfc = J drV^^[d^duip*ipk]- Noting that the 
orbitals ipj form an orthonormal set we see that {da — T]a)'4^k — Yli ^ajki^j with similar 
relationships for the others. 

This allows us to write the matrix eigen-equation 

Q^NPiv =RPjv (91) 
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where P^v is an A^Bas x 1 matrix with elements pnj- Here 




Orthogonahty is given by PJ^NFm = (^atm- We note that all tensor components are bundled 
into N and M via summation over a, fi and u. These equations are only true in the strict 
limit of an infinite number of orbitals but will converge to the true approximate solution in 
typical systems. 



C. One-electron like asymptotic behaviour 

In the asymptotic |r| — > oo region of a localised system, the highest occupied molecular 
orbital (HOMO) ipH with occupation number fn dominates [at least in the non- degenerate 
case or for degenerate system with spherical symmetry where wc set fn = {21 h + 1) and 
i^nix) = -Rnir(^)) the radial component of '0ji^/^j„(r) = i?„^(r)y;^^(f2r)]- Thus, in the limit 
r — >■ oo we can set: 

^fH\M\ V (94) 

t^,u^-\d^ri,, S^2{V'^^-eH). (95) 

All equations reduce to the one-electron form in these limits, with the single occupied orbital 
given by the HOMO. 

We can check how close a region of a many-clcctron system is to a onc-electron-like system 
by considering the difference of S from its asymptotic form. Here we set 

S{r) =S{r) - 2{V'^^{r) - e^) (96a) 
=C°(r) + 2(e^-e») (96b) 



where 



^V)^E^^51f^-i.wr^ (97) 
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and where C° is related to the Fermi- hole curvature^^ of the KS system. In the limit |r| — > oo 
it is clear that C°(r) = and e(r) = €h (true everywhere for one-electron systems). 

This outer region is crucial for van der Waals physics in many systems. Typically only 
transitions within a small frequency range, especially that to the lowest unoccupied molecular 
orbital (LUMO), will dominate in this region. In such transitions J driplSipk will be small. 
The asymptotic 'exactness' in this limit suggest strongly that the scalar approximation to 
CM is a appropriate for calculating vdW forces. 

Even a well chosen but simple approximation, using a limited number of unoccupied 
orbitals and quasi-orbitals, may provide quite accurate estimations of the vdW physics. 
The ability of the CM (and the scalar approcimation) to make collective modes for the 
transition physics may aid in convergence compared to full dRPA calculations. 



D. Transitions and energy calculations in the scalar approximation 

Once we have obtained solutions of (91) we can use pNj to evaluate djv through (80). 
Using also we can then evaluate Xo and properties which depend on it such as the 
correlation energy. 

We can also use the asymptotic form of S{r) to further investigate and through it 
correlations. Combining (80) and (96) lets us write 



=*°^p^, (2^H + ^)^, (99) 
j 

where the fljH — €j — en comes from applying the Schrodinger Hamiltonian —\V'^-\-V^^ to 
. Here S covers the deviation of many-electron systems from their one-electron equivalents. 
Let us expand (99) as follows: 

dN ='^^^dNj^j. (100) 
3 

We define matrices O with Ojk — Sjk^kH and § with Sjk — J dr'0*'§V'fe- Thus d^j — [^N]j 
where Dat = ^20 -|- Pat and Pat is a solution of (91). We can use this expansion to 
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calculate the Coulomb projections Wnm defined in (37). Using (100) we find 



Wnm ^B^W^^^Bm 

=p1^(20 + §)W('^)(20 + §)P7v (101) 

where 

We can use Wnm to calculate the correlation energy through equation (38). 

In molecular orbital language it is clear the calculation of W^'''^ is the only step involving 
four-center integrals required by a correlation energy calculation. However, in contrast to the 
four orbital (two occupied, two unoccupied) terms {ia\jb) required by a full KS calculation it 
requires two uncontracted orbital indices (occupied or unuccopied) only. For large molecules 
this is a substantial saving. 



IX. ID HARMONIC OSCILLATOR 

We define a ID non- interacting and spin-free Harmanoic oscillator with N^, electrons in 
the groundstate and V^^{x) = x'^/2. Here the KS orbital wavefunctions and energies take 
the form 

^j{x) ^K,Hj{x)e-^'/', ej (103) 

where Hj{x) are Hermite polynomials with j > and Kj — [-^/7r2-'j!]~^/^. Since the orbitals 
are filled up to j < A^e — 1 the density is 

and it is clear from the expansion of H^^-i that Iim,.^oo 'n^{x) ~ r^'i(N^^~i)^-x^ lim2;^.oo vi^) 
—X. Thus the system is well-bounded and the KS-orbital approach should work well. 

In the one-electron case 77 = — x, t — \ and V^xx — 1- and we can solve everything 
analytically. The matrices required to form N and IR in equations (92) and (93) have elements 
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Tjk — |(^jfe, Vjk — ^jk a-nd 

Djk =i^ji^k J dxe-^"/^Hj{x)[{d, + x)Hk{x)e-^"/^] 
—2kKjKk J dxe"^^ Hj{x)Hk_i{x) 

Here [DtD],-^ = 2jSjk, P^B>%k = - l)5,k and p^^J^%k = - l)(j - 2)5jk. 

Thus Njk = 2jSjk and Rjk = 2jdjk + 6j(j - 1)5,^ + 2j{j - - 2)djk = 2fSjk. We 
can ignore the j — solution as it will not contribute to xo- We thus choose solutions with 
Nj > where p^.j — b]<j.j 1^/22. Here — — | = j as expected and the transition 
density mode is dfq.{x) = -ilJo{x)[{d'^ + 1 - x'^)'4)j{x)]/y/2j = y/2nM,i^oix)iJj{x). 

In the many-electron case we can solve the problem semi-analytically for much of it, 
requiring numerics only for terms involving and subsequent diagonalisations. Errors 

occur due to truncation of the basis set but integrals can be obtained with near exact 
accuracy with Gauss-Hermite quadrature. It should be noted that the two-electron system 
is not predicted exactly by the CM as there is no spin degeneracy. 

The rapid convergence of the method is demonstrated in Table I where we show the 
convergence of the fourth transition frequency ^4 for the two-electron system. With as few 
as 15 states the error is already under one part in ten thousand. 
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3.8801 


3.8802 


3.8802 


logio |Err| 


-1.22 


-3.95 


-4.05 


-5.06 


-7.49 


—00 



TABLE I. Convergence of the fourth transition frequency for the two electron system. Here Err = 



The KS transition frequencies of an A^e electron system have energies Qai — — = J 
and each frequency J has multiple contributing transitions. Here min( J, N^) modes have 
transition frequency J and the density transitions are proportional to ipi{x)ipi+j{x) when 
i < iVe and i J > A^e- 

We present the CM transition frequencies in Table II for systems with up to 20 electrons. 
The CM has single-valued frequencies distributed approximately with the integers. Each 
CM transition density is therefore composed of multiple KS-like transitions which, by the 
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sum rules (32) must have their weights haiN dominated by transitions with D,ai ~ JIat- In 
Table III we show some weights h^N for the five electron system. It is clear that the A'"th 
CM transition puts most weight on the HOMO-HOMO+A^ [5 - (5 + N)] KS transition as 
one would hope. 
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TABLE II. First eight, 10th and 20th distinct eigcnfrcquencies of the KS system and CM system 
with different numbers of electrons N^,. The KS transitions are all integer and degenerate with 
min(il, A^e) transitions for a given integer frequency. 



It is remarkable that the lowest three frequencies are exact to four decimal places for 
iVe > 1, particularly as the second smallest is comprised of two transitions and the third of 
three in all but the two-electron case. By the sum rules (32) this means that each mode must 
be comprised solely of transitions with the given frequency, and as a consequence, both VLjq 
and (In are predicted exactly for the HOMO-LUMO transition. That such a relationship 
holds for as many as 20 electrons demonstrates a strength of the physics and a resihence of 
the approximations in the CM approach. 
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X. CONCLUSION 



In this paper we have reformulated, simphfied and investigated the CM formahsm of 
refs 4-7, and provided a direct proof of its exactness in one-electron systems. We have 
described two formulations: a normalised approach very similar to the usual form, and a 
scalar approximation that is directly equivalent to the usual tensor form for one-electron 
and one-dimensional systems. 

The orthonormal form described in Section VI is vital for the description and calculation 
of localised systems, where the standard w-based formulation does not behave well. This 
reformulation is required for calculation of atomic, and molecular systems using standard 
basis set approaches involving GTO and STOs. It also provides a relatively simple way of 
proving (see Section VII), direct from the Schrodinger equation, that the CM is exact for 
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TABLE III. Weights haiN in % with > 1% contribution for the five electron system. The KS 
transition frequencies are Q^ai = a — i. Tabulated weights may not sum to 100% due to the absence 
of < 1% contributions. 
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one-electron systems. 

The scalar approximation of Section VIII provides a simplified approach to the CM, with 
some sacrifice of some exact properties of the full theory in general systems. It could be 
used as a faster alternative itself or as a doorway to further and cruder approximations. In 
particular it opens up analysis of transition frequencies and densities in a simple manner 
by relating information from the KS orbitals to the collective modes in the CM. The ap- 
proximation is exact in one-electron systems and should accurately predict the physics of 
asymptotic regions where the behaviour is essentially one-electron like. 

In this approximation, at least for the difficult four-center term required by correlation 
energy calculations, the scalar CM approximation has an efficiency gain of 0{Nf) over 
regular KS where Ne is the number of occupied electronic orbitals. This means that in 
large molecules (and in periodic bulks), evaluation of Ec under the approximate CM will 
be a significantly quicker calculation, with corresponding lower memory requirements. The 
approximation also makes a gain of 0(3^) in diagonalisation compared to the full CM theory 
and 0(3^) in the four-center terms. 

Both reformulations presented here will aid in the development of a basis set based, CM 
approach for atomic and molecular systems. The difficulty in such an approach is the pres- 
ence of non-constant denominators in t^,^, rja and S. We believe, however, that a tractable 
approach to dealing with these should be possible by using Gauss-Hermite quadrature with 
GTOs. 

It should also be possible to further develop vdW functionals that use a further approxi- 
mation to the CM based on cleverly chosen orbital-like basis functions, possibly incorporating 
some unoccupied KS orbitals. Such functionals could incorporate ideas from refs 15-17 to be 
made even more efficient than the functional described here or in ref. 7, without sacrificing 
the vital physics that make the CM attractive. 
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Appendix A: Kinetic stress tensor 

There is some debate over the appropriate form of the kinetic stress tensor T. We have 
chosen the form given in (3) and (4), whereas the original CM papers used the form given 
in (13). One requirement is that, for Kohn-Sham systems, T° obeys the force balance 
condition V • T° = —nPS/V^^ which is true both for the form (4) used in this manuscript 
and for equation 17 of ref. 5 [leading to our equation (14) for a Kohn-Sham groundstate] . 

There are various ways to arrive at the possible forms of T(r) (see e.g. refs 29-33). 
For completeness, we show one way to arrive at the form (3), (4) chosen for the present 
work. This form helps to make the CM equations particularly simple. We base our form on 
the Wigner transform of the "classical" kinetic stress tensor. In a "classical" fluid with no 
current present, T is defined by 



Here / is the classical one-body distribution function - i.e. the phase space probability 
density for finding a particle at position r with momentum p. To model a quantum system 
we can replace f{r,p) by the Wigner form 



where p{r, r') is the quantal one-body density matrix. Using the short-hand p = p{r -\- 
x/2,r — x/2) we thus find 




(Al) 




(A2) 




(A3) 



and so, by the chain rule. 




dx^d^,p{,r + x/2,r - x/2)\ 



x=0 




(A4) 
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with df^ = and = dr'^. Furthermore 

= ^5M^i.p(T','r) = ^9^9^n(r) (A5) 
Combining (A4) and (A5) we find 

T,.{r) =\ [{d,dl + d'^dMr.r%^^, 

- \d^dMr)- (A6) 

for a state with zero current, in agreement with (3). 

Since p{r,r') = Xlt ^ groundstate KS system we find the stress tensor 
T° for this situation to be 

i 

- \[d,d.n'{r)] 

=3^E/^[^^^*H9.V'«] - \[d,d.n'] (A7) 

i 

Thus (3) is exactly the Wigner version of the "classical" kinetic stress tensor (Al). 
Using (5) we demonstrate that 9qT°^ = —nPV^^. Here we note that 

9aTl = lj^f,tj, (A8) 
j 

where 

and we have used the Schrodinger equation V^ipj = 2(y^^ — ej)ipj to derive the final ex- 
pression. Finally d^T^a — ~ Ylij /j^il^^^lV'jP — ^'^^ the proof is complete. 

Comparing this with the form T° (14) used in the original CM formulation, we find 
daT^ — daT^a siuce daid^danP] — [9^V^n°] = d^Sij^lV'^nP]. Since these are precisely the 
terms that differ between (4) and (14) it follows that the gradients must be identical. 
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Appendix B: Compact form of the CM equations 



In an earlier work^ we state without proof that equation (8) defined via (11) here (equa- 
tions 2-5 of ref. 7) are equivalent to equations 14-16 of ref. 5. The only non-notational 
difference is in the kinetic force term ir^^'" — — Kn defined here via (9), in equation 14 of 
ref. 5, and in equation 53 of ref. 6 (abbreviated as G53) where it is derived from their 
equation C8 (abbreviated as GC8). 

We demonstrate that the two forms are equivalent by working from GC8. While the 
same result can be obtained directly from G53 the derivation is less clear and less succinct. 
Specifically we must show that 

=^--K,^u^ (Bl) 

where is defined in (9), since the remaining terms in (8) follow directly from equations 
14-16 of ref. 5. 

Following Gao et al^ we write GC8 as 72 [w] = / dri where 

+ ^f^/xi/^^C^aa + -^Uo,,^dyU^^aj, (B2) 

with U^y = l{Uf,^u + Uu,^,), Uaa = Ua,a, and ]Cf,u = I'^Y^i fi[d^^i^*][duipi]. 

We can use integration by parts to remove derivatives of n° from the integrand I. As 
such 



If/ 



dr 



71° 



(B3) 



where we define 7"°^ = S/C^ui, — 1^1°^^^ to be equal to (4) rather than the form appearing in 
refs 5 and 6. The undefined terms in the integrand of (B3) take the form 
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where we have expanded U/^i, and used the product rule on all derivatives to arrive at these 
forms. The following identities are also used in the derivation of Z: 

which follow from exchange of indices under summation (eg. Aj^fj^B^i, = A^j^Bi^ij). 
We can expand the terms of y^i^^ni, as follows: 

where we again exchange indices where appropriate. This leads to the following result 

and thus yixv,ij.v — Z — Uu,fj,aUij,,i^a- The cancellation of so many terms is quite remarkable. 
Finally (B3) becomes 

} (B8) 
j (^ru^K^^u„ (B9) 

+ daTj^^da + d^Tld^ + d^Tl^d^. (BIO) 
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Here we used integration by parts on the derivatives of to obtain (B9) and (BIO). The op- 
erator K defined in (BIO) is identical to that defined in (9). Taking the functional derivative 
w.r.t. U/^ thus gives 

=^ - -^M.^. (Bll) 

and it is clear that (Bl) is satisfied. 



Appendix C: Correlation energy expressions 



Let us first define the projected Coulomb operator Wnm = {um\Q\un) where Q is 
defined in (16) or (17). This can be written as 

J drul^ir) ■ [QunKt) (CI) 

= J drdr'<(r) • Q(r, r') ■ u^ir') (C2) 



I 



c\vc\v' 

p-^dl,(r)d^(rO. (C3) 

where similar equivalences hold true for Wnmx — {umx \ Q \un\) or for alternative forms of 
the Coulomb potential (e.g. range-separated). 

Using equations (18)-(19) and (24)-(25) we can write the correlation energy as 

.J. TdA r^^ ^d-"-' 



2 Jo Jo 



^ ^ Wnnx Wnn 



2 .In \ ^VInx 20iv 



(C6) 



Here the governing eigen-equations for unx are as defined in (19) 

nl^UNX =(R + XQ)UNX (C7) 

with normalisation J 'nP{r)u*j^^{r) • UMxif) — ^nm- 

Following the ideas of Furche,^^ we can take the A derivative of J drwjvA-(C7) to work 
directly from (C6). Here 

2^Nx[dx^Nx] =jKx- Q^iVA- (C8) 
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where derivatives of Un\ can be ignored by the Hellman-Feynman theorem. Thus 

{r)[QuNx]{r) _ Wnnx 



[dx^Nx] = J dr- 



2Qnx 

and we can write the correlation as a sum over zero-point energies such that 



N 



Oat-^^jv 1 + 



NN 



2n% 



(C9) 

(CIO) 
(Cll) 



where Q,n — ^m- In certain systems it may make sense to work with exchange and corre- 
lation together. Here 



^xc =1 Yl ~ ~ ^ / dr-n°(r)«;c(r) 



(C12) 



where wc{r) — J dr' ^^^_^,} is like the Coulomb potential at zero distance. 



To solve for involves a difficult diagonahsation of R-l-Q and may be best avoided. We 
can use the complete and orthogonal nature of Unx and Un to write Unx = '^k^nkUk 
where U^U = I. Thus we can write 



Xx{r,r') = - ^'YXNMxd*N{r)dM{r') 



(C13) 



NM 



where X^mo = Snm/{o''^ + ^n) or Xq = (cr^ + L) ^ where Lnm = ^nm^n- Solving for 
Xa = Xo + Axo^Xa gives 

I 



This corresponds to 



X, 



(72 + L + AW 



(C14) 



X Tr 



W 



W 



(C15) 



_(72 + L + AW (72 + L_ 
which can sometimes prove useful in real calculations 

We can relate all this back to (Cll). Solving the eigcn-cquation VD^ = (L + AW)V gives 



^MiVA=[V((7'+BA)-V]MiV 



(C16) 

(C17) 



K 
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With appropriate sorting of the eigen- value/ vector pairs it is clear that the eigenvalues Dattvi 
of L + W are Qj^ — Q^f^ and U = Thus the two expressions are equivalent. 

Furthermore for » Wn where Wn = I^m I^a^mI we can solve the eigen-problem 

pcrturbatively so that Q%^n% + Wnn- Thus = \/^^TW^ ^^1^(1 + ^ - 

We can speed calculation and improve numerical stability by choosing an N* above which 

we use the approximation. Setting — ^^r- we then find 

J2 [^N - lliv (1 + M] - E (C18) 

N=0 N>N* 

Such a perturbative approach will be almost guaranteed convergent if we replace the 
Coulomb potential by the range-separated^*^"^^ form Vq^^\r) = erf(gRs-R)/-R. 

Appendix D: One-Electron governing operator 

As discussed in section VII [equations (61)-(67)] we can obtain the exact hnear pertur- 
bation solutions in a one-electron system via a solution of 

u'^ =R^^"£ (Dl) 

where the associated change in density is n 1 = - WC^' . Since this is an identical problem 

to finding the CM solutions via 

=R«)£ (D2) 

it is clear that if R^^ = R^^^ then the CM is exact for one-electron systems. 

For a one-electron system we set V = V^^ — such that the groundstate occupied orbital 
is a solution of [— |V^ + V]'?/' = 0. Equations (55) thus become rja = <9Qlog IV^I and S = 2V 
and thus h — — |(V^ — 5"). The former can be used to derive the following 

where these identities will be used throughout this appendix. 

Comparing the final term of (67) with (11) we wish to swap [3^ — 7]^) and [8^ + 7]^) across 
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sides. This can be done by repeatedly using 



dj =fd, + [dj] 
vV=/v' + [vV] + 2[a„/]a„ 

=/V^-[VV]+29„[9«/] 
fV'=V'f + [V'f]-2d^[dJ-] 

=vV-[vV]-2[a,/]9, 



in the appropriate places. 

We first address the terms involving 5". These are 



and thus 



(D3) 



(D4) 



The terms involving the Laplacian are a httle more difficult to deal with. Here 
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and 



- 2ri^[darin]da + 2da[dariu]rin 
Combining the above lets us write 

- 4(9,. + r]„)tan{da - rja) 
+ 2[a„7/,][a„r7^] -2[VV] 

We can now use (D6) and (D4) together to show 

A ple~ E>le~ pie" 

- 4(a„ + ria)t^,u{da - Va) 

- 4(9^ + r]^)tat,{da - rja) 
+ S^^^ + 2[dar]u] [daVt^ 
-2[VV] -4[?7a9«V] 
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where we cancel most terms via 



+ [VadJuu] + St 

[^uV^^S]^-2St^, + [r)^S,,]. 

The operator terms are now the same as those of (50). The remaining constant is 

4^K^. ^S,^, + 2[daV.][daVn] - 2[V V] - "^iVadJ^.,] 
- 2[d^r]^d^T]^] + 2[r]adadf,'qy\ 

where we have used [V^i^j.] = -^S^^i, + [dt,rja][d^rja] + [r)adi,dij,r)a\ and [daVu] = [di,r)a] to 
arrive at the final form. 

Thus (D8) becomes 

- {da + r]a)tiuu{da - Tja) 

- {da + 'na)tau{d^ - T]^) 

- {du + Vu)tan{da - Va) (D9) 

which is identical to R^^J and through it R^iy via (45). 

Finally, since t/j^ and through it 0^ and n} can be calcuted from ^ it is clear that a solution 
to the CM is identical to a hnearised solution of the Schrodinger equation for a one-electron 
systen, and vice versa. 
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Appendix E: Scalar CM operator in the KS basis 



Inserting (50) [and setting ( — S) — {da + rja) {da — rja)] into (85) and using integration 
by parts we find 



It follows from rj = ^Vlogn^ that [d^r]i,] = [dur]fj] and thus dfj,r]u + r]fj,du = 'qi.d^ + di^r]^. 
Therefore D^Dv = D,,!)^ and we can simplify (El) (noting that we can also swap Greek 
indices as they are summed over) to 
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